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Abstract. The asymptotic tails of the probability distributions of thermody- 
namic quantities convey important information about the physics of nanoscopic 
systems driven out of equilibrium. We apply a recently proposed method to an- 
alytically determine the asymptotics of work distributions in Langcvin systems 
to an one-dimensional model of single-molecule force spectroscopy. The results 
are in excellent agreement with numerical simulations, even in the centre of the 
distributions. We compare our findings with a recent proposal for an universal 
form of the asymptotics of work distributions in single-molecule experiments. 
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1. Introduction 

Traditionally, statistical mechanics is concerned with averages; the probability 
distributions for thermodynamic quantities of macroscopic systems are so exceedingly 
sharp that only their most probable values matter. These in turn are practically 
identical with the averages. Only near instabilities, deviations from averages become 
important. In static cases these fluctuations are well described by the second moments 
of the respective distributions; if the dynamics is of importance as well, they are 
characterized by the correlation functions. 

When investigating nanoscopic systems from the point of view of thermodynam- 
ics, the situation changes. Fluctuations are now strong and ubiquitous, and corre- 
spondingly, the probability distributions of relevant quantities are broad and poorly 
characterized by their leading moments alone. Whereas this seems rather obvious, it 
came as a real surprise that work [T] and fluctuation [21 [3] theorems which form the cor- 
nerstones of the emerging field of stochastic thermodynamics [4j [5j El O [8] are relations 
that probe the very far tails of the respective probability distributions. Very unlikely 
events now carry significant information about the physics of the system under consid- 
eration. On the one hand, interest in mathematical investigations like large deviation 
theory (see [9] and references therein) is renewed, on the other hand, techniques are 
in demand that allow to determine the asymptotics of probability distributions. Since 
rare events are hard to get in experiments and numerical simulations, approximate 
analytical procedures have to be developed. 

In the present paper we apply a recently proposed method for the analytical 
determination of the asymptotics of work distributions in driven Langcvin systems 
|10j to a simple model for single-molecule force spectroscopy. In these experiments 
(for a recent review see [11]) the free-energy difference A.F between the folded and 
the unfolded state of a biomolecule is determined from the distribution of work W 
obtained in isothermal unfolding and refolding processes. If only one transition is 
monitored |13j . the Jarzynski equation [I] 

( e -^) = e ^ AF (1) 

is employed, where (3 denotes the inverse of the temperature. If histograms of work for 
both the forward and the reverse process are compiled |14j . it may be advantageous 
to use the Crooks fluctuation theorem [15] 



Pr(~W)- e ■ 

Here, the free-energy difference AF is identified from the intersection of the probability 
density functions P(W) of the forward and P r (—W) of the reverse process. Note that 
in both cases an accurate estimate for AF requires reliable information about the tails 
of the work probability distributions. 

The paper is organized as follows. In section 2 we introduce the notation 
and review the general method |10j . Section 3 contains the application to an one- 
dimensional stochastic process with a time-dependent double-well potential and the 
comparison with numerical simulations of the system. In section 4 we discuss our 
results and in particular compare them with a recent proposal for the general shape 
of the tail of the work distribution in single-molecule experiments |19j . 
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2. General Theory 



To pave the way for the analysis and to fix the notation we summarize in the present 
section in a very condensed way the main steps of our approach to determine the 
asymptotic tail of work distributions in driven Langevin systems. For more details 
the reader is referred to [10] . 

We investigate a system described by an one-dimensional, overdamped Langevin 
equation which in dimensionless units has the form 

x = -V'(x,t) + y/2/PZ<fi- (3) 
Here x denotes the degree of freedom, V is a time-dependent potential modelling the 
external driving, and is Gaussian white noise with (£(t)) = and = 
S(t — t'). We denote derivatives with respect to x by a prime, and those with respect 
to t by a dot. The initial state x(t = 0) =: xo of the process is sampled from the 
equilibrium distribution at t = with inverse temperature /3, 



Po(xo) = — exp[-pV (x )] 



(4) 



(5) 



Accordingly, Vo(x) := V(x, t = 0) is the initial potential and 

Zq = J dx exp [— (3Vo(x)} 

the corresponding partition function. 

The work performed on the system for a particular trajectory x(-) is given by |16j 

W[x{-)} = f dtV(x(t),t) . (6) 
Jo 

Due to the randomness inherent in x(-), the work W is itself a random quantity and 
its pdf can be written as 

P(W) = J^Tr exp[-pVo(x )} J dx T 

x(T)=x T 

Jvx(-) p[x(-)} 5{W-W[x{-)]) 

x(Q)=x 

For mid-point discretization we have |17j 

rT 



(7) 



P[ x (')] = •^[ :^; (•)] ex P 
with the normalization factor 
Af[x(-)] = exp 

Hence 



-~ [ dt (x + V'(x,t)y 
4 Jo 



- / dtV"(x(t),t) 



(8) 



(9) 



dq 

4tt//3 



x(T)=x T 

r Vx(-)N[x(-)]exp{~(3S[x(-),q}} 

x(0)=xo 



(10) 



Asymptotic work distributions in driven bistable systems 



4 



with the action 

S[x(-),q\ = V (x ) + f dt 
Jo 



Y \w. (11) 



We are interested in the asymptotic behaviour of P(W). Rare values of W are brought 
about by unlikely trajectories x(-). In the spirit of the contraction principle of large 
deviation theory [9], we expect that in the asymptotic tails of P(W) there is one 
trajectory for each value of W that dominates P(W). To find it, we have to maximize 
P[x(-)\ := po(xo)p[x(-)] under the constraint W = W[x{-)\. This can be done by using 
a saddle-point approximation of the integrals in (fl"U|) . The result is 

P W = ^!^tfi(l + 0(l/«). (12) 
V R dot A 

To use this expression in explicit examples, we first have to determine the most 
probable trajectory x(-) satisfying the Eulcr-Lagrangc equation (ELE) 

S + (1 - iq)V - V'V" = (13) 

where V(t) := V(x(t),t) and similarly for derivatives of V. The ELE is completed by 
the boundary conditions 

i - Vq' = 0, S T + V T = (14) 

and by the corresponding value q of the Lagrange parameter q ensuring W[x(-)] = W. 
Using x(-) and q, we calculate S := S[x (■),(}] and N := J\f[x(-)}. Then all terms in 
(|12[) depending solely on the optimal trajectory are determined. 

The denominator Vi? det A in ([T^|) comprises the contribution from the 
neighbourhood of the optimal path and stems from the integral over the Gaussian 
fluctuations around x(-) and q. Here, 

A := + {V") 2 + V'V'" - (1 - iq)v" (15) 

at z 

denotes the operator of quadratic fluctuations which acts on functions <p(t) on the 
interval < t < T obeying the boundary conditions 

v^o) - m - o , 

Vg<p(T) + <p(T) = . 

A simple prescription to calculate det A is as follows [18] . Solve the initial value 
problem 

A X (t) = , 



(17) 

X(0) = 1, x(0) = ^ " 



then 



detA = 2(V^ X (T) + x(T)) . (18) 

The factor R in (TT^I) accounts for the influence of the constraint © on the 
fluctuations around x(-). Since also the trajectories from the neighbourhood of x(-) 
have to yield the very same value of W, fluctuations violating this constraint are 
suppressed. This gives rise to a correction to the fluctuation determinant of the form 

R= dtV (t) A- l V (t) , (19) 



Asymptotic work distributions in driven bistable systems 



5 




Figure 1: Evolution of the potential V(x, t) ([22]) in the time-interval < t < T for two 
exemplary parameter sets. Shown is V(x, t) for t = 0, t = T/2 and t = T respectively, 
(a) a = 0.25, b = 0.5, c = 2, r = 3, T = 6 and (b) a = 0.025, 6 = 0.5, c = 3, r = 1, 

r = 9. 



where ^4 1 denotes the inverse operator of A. To explicitly determine R, it is 
convenient to solve the ordinary differential equation 

Aip(t)=V'(x(t),t) (20) 

with boundary conditions (|16| and to use 

R = / dt i/j(t)V (t) . (21) 



3. The driven double- well 

The unfolding and refolding of single molecules can be modeled by a time-dependent 
double- well potential of the form 

V{x,t) = ax 4 - bx 2 + r(c ~ t)x , (22) 

where x denotes the extension of the molecule in the direction of the force p~3l [Ml [T9] . 
The parameters a and b characterize depth and separation of the two minima of V, 
c fixes the moment at which V is symmetric, V(x) = V(— x), and r denotes the 
transition rate. Choosing T > c for the final time, the two minima will interchange 
global stability during the process. We have used two exemplary sets of parameters 
for which the time evolution of V(x,t) is sketched in figure [TJ The main differences 
between the two sets are that in (b) the minima are further apart and the transition 
rate r is smaller. 

We now apply the analytic method to obtain the asymptotics (fT2"j) of the work 
distribution of the dynamics defined by the potential V(x,t) (f2"2"]) . To clarify the 
procedure, we compile the necessary equations. We have from (JS|) and ^ 

Z = J dx cxp [~/3(ax 4 - bx 2 + rex)} (23) 

M[x{-)] = exp 6a f dt x 2 (t) - bT . (24) 
Jo 
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Figure 2: Trajectories for the potential shown in figure [Ta] for (a) the forward and (b) 
the reverse process. Shown are optimal trajectories x(-) for exemplary work values 
W . For a comparison is plotted on top (thick lines), the average trajectory of the 



simulation (full blue line), and the average 



\ eq 



from the equilibrium distribution 



corresponding to the instantaneous values of the parameters (dashed line) . 



The ELE ([13]) reads 
,2-5 

— OZUX" U -I- 71 C — U I I iiUI — ZUI 

(25) 



- 32a.T 3 6 + r(c - t) [\2ax 2 - 2b] 



+4b 2 x + r(l - iq) , 

and its boundary conditions (|14[) are of the form 

= xn — 4axn + 2bxn — rc , 

° (26) 
= x T + <iax T - 2bx T + r(c - T) . 

The constraint © is 

W = -r [ dtx{t;q) . (27) 



Using a standard relaxation algorithm, the numerically solution of equations (|25[) - 
(|27| for desired work values W yields optimal trajectories x(t) depicted exemplarily 
in figure [21 

The operator A from (|15p acquires the form 
d 2 

A = - — T + 240a 2 x 4 - 96abx 2 

at' 1 (28) 

+24ar(c - t)x + 46 2 
with the boundary conditions (|16|) 

0= [I2ax 2 -2b]^0)~m, 
= [\2ax 2 T - 26] ip(T) + <p{T) . 



To obtain det A according to ([18)) , we determine 

det A = [2Aax?r - 46] X (T) + 2x{T) (30) 
by solving numerically the initial value problem (|17p 
X(t) = [240a 2 x 4 - 96abx 2 

+24ar{c - t)x + 4b 2 ] X (t) = , (31) 
X(0) = 1 , x(0) = I2ax 2 - 2b . 
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Figure 3: Trajectories for the potential shown in figure [Ta] for (a) the forward and 
(b) the reverse process. Shown is, exemplified by single realizations &(•), the range of 
trajectories attainable from the simulation (full grey lines), the average trajectory of 
the simulation (full blue line), and the average (x)° q from the equilibrium distribution 
corresponding to the instantaneous values of the parameters (dashed line) . 



This has to be done for each value of W separately by using the appropriate results 
for x(t;W) and q(W). 

The last ingredient for the pre-exponential factor is R from ([21]) . To determine 
it, we need to solve the boundary value problem (|20|). (fl6|) 

i'(t) = [240a 2 x 4 - 96abx 2 

+24ar(c - t)x + 4b 2 ] ip(t) + r , 

= [l2ax 2 - 26]V(0) - ^(0) , 

= [\2ax% - 2b] ip(T) + ip(T) 

for each x(t; W) and q(W) and use the result in (|2"Tj) 

R=-r [ dt ij){t) . (33) 
Jo 

Plugging the numerical results for Z , TV", S, detA and R into (|12p . we obtain 
the final result for the asymptotic form of the work distribution. We carried out this 
program for the two processes characterized by the parameter sets given in the caption 
of figure [TJ including for both cases the reversed processes defined by the substitution 
t — > (T — t). From simulations of the Langevin equation ([3]), we also obtained from 
([6]) the corresponding histograms of work distributions. The values of j3 and T are 
chosen such that most of the trajectories reach the right minimum of V at the end of 
the forward process, i.e. the molecules are stretched until virtually all of them unfold 
(cf. figure [3]). The results for the asymptotics are shown in figure HI together with the 
outcome of our numerical simulations. 

In a recent paper |19j , Palassini and Ritort propose an universal form for the tails 
of work distributions for single molecule stretching experiments given by 
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Tabic 1: Fit results for (|34[) for the two examples (a) and (b) shown in figure [T] for 
the forward (fw) and reverse (rv) process. TV denotes the number of realizations used 
in the simulation, the other parameters are from the proposal (|34l) . 



Set N n SI a W c 8 A 



(a) fw 10 4 4.27E-4 9.31 -26.4 11.5 3.35 4.22 

(a) rv 10 4 1.16E-5 8.64 -31.2 38.2 3.26 2.26 

(b) fw 10 5 4.97 8.61 -3.01 19.9 2.47 1.45 
(b) rv 10 5 1740 26.1 26.5 61.2 3.33 0.607 




Figure 4: Work distributions P(W) (forward process) and P r (—W) (reverse process) 
for the two potentials (a) and (b) shown in figure [TJ Corresponding results obtained 
from numerical simulations are shown by open circles (forward process) and open 
squares (reverse process). Full lines depict the asymptotics according to (fl"2"j) . dashed 
lines are fits of (f3"4"|). The value of the free-energy difference AF obtained from 
numerical integration of the partition functions Zq, Zt (|23p is indicated by the vertical 
line. 



Here, W c is a characteristic work value, 51 > measures the tail width, and n, a 
and 8 are further constants. The Jarzynski equation (TTJ) stipulates S > 1. Based on 
(|34[) . Palassini and Ritort present in [TH] three slightly different analytical methods 
to improve the estimation of free-energy differences AF from the Jarzynski equation 
(QJ. To decide which approach to use, they distinguish three regimes defined by the 
parameter combination 

A := ((5/0) 5/( ' 5 ~ 1) IniV . (35) 

The three regimes are then identified by A > 1, A <C 1 and A < 1 respectively. They 
test their method with experimental data from DNA stretching experiments. For more 
details regarding the improved estimation of AF and the experiments see |19j . 

We fitted the empirical asymptotic form (|3"4")l to the tails of the work distributions 
obtained in our simulations by standard least-square fits, starting with a Gaussian 
distribution specified by W c = (W), tt 2 = 2 ({W - W c ) 2 ), a = 0, 6 = 2 and n = l/^/n. 
A subtle point in the procedure is to find the optimal interval of work values for the 
fit. The resulting parameter values are listed in table [TJ the corresponding fits are 
included into figure 01 
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4. Discussion 

As shown in Figure H] for both parameter sets (a) and (b) , we achieve an excellent 
agreement between simulation and asymptotics, not only asymptotically for the tails, 
but also for the centre of the work distributions. The good reproduction of the centre 
of the distribution is presumably due to the fact that also the probabilities of typical 
work values are dominated by single trajectories and their neighbourhood. Note also 
the perfect match between the free-energy difference AF and the intersection of our 
asymptotics PiW) and P r {— W), which demonstrates that the Crooks relation @ 
holds in its exact form for our asymptotics, as has been shown in [TO] . 

In addition, figure [5] illustrates the optimal trajectories x(-) which dominate 
the asymptotics of the work distributions (|12[) . In comparison with the trajectories 
obtained from our simulations shown in figure [3j the trajectories i(-) contributing 
to the tails of the work distribution arc of much broader variety. Interestingly, the 
probability of small work values is dominated by x(-) that run into the evolving right 
minimum even before the minimum is shaped. 

For the parameter set (a), the fit of (|34|) compares well with the analytic 
asymptotics. In this case a combination of a histogram from experimental values 
and (f3"4"]l would therefore result in reliable estimates for the free-energy difference AF. 
For the parameter set (b) only the forward process is well described by the fit, for 
the reverse one the tail of the distribution is markedly overestimated. This shifts the 
intersection point between P(W) and P r {— W) away from the correct value of AF 
as can be seen in figure I4bl Also, some fit parameters for this case listed in table 
[TJ clearly deviate from the other cases. To investigate that mismatch more closely, 
we also fitted (|34|) to the analytic asymptotics (not shown). This results into similar 
conspicuous fit parameters, but the fitting curve now is almost congruent with the 
analytic asymptotics. From that we conclude that the empirical asymptotics (|34[) is 
also valid in this case, but the number of work values obtained from the simulation 
is not enough to reliably fit ([34]) . Note that rather than taking the intersection point 
between P(W) and P r (—W), Palassini and Ritort use in [19] much more sophisticated 
methods to obtain estimates of AF, based on their empirical asymptotics (|34| . Our 
investigation aims only at validating (|34l) . 

In addition to the two parameter sets displayed in figure (Q}, we also investigated 
several other realizations of the potential ([22]) (not shown). Mostly, we found that fits 
of (]34| to histograms of 10 4 to 10 5 work values extrapolate reliably to the asymptotic 
tail of the distribution. But as exemplified by the case shown in figure 14b] there is no 
guaranty that a number of 10 5 work values is sufficient. 

5. Conclusions 

The asymptotics of work distributions for driven Langevin systems can be determined 
by the method proposed in [10]. We employed this method to determine the 
asymptotics of a system defined by the potential (|2"2"]) modelling the stretching of 
single molecules. The unfolding of the molecules corresponds to the forward process, 
the refolding to the reverse process. We obtained histograms of work by simulating 
the Langevin equation ([3]). The form of the work distributions are found to be near- 
Gaussian, similar to distributions measured in a DNA stretching experiment [19] . We 
observe excellent agreement between asymptotics and work distribution, not only for 
the asymptotic regime but also for the whole range of work values. 
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One aim of single molecule experiments is to obtain the free-energy difference 
between the folded and unfolded state of the molecule. If both the work distribution 
for the forward and reverse process is available, the Crooks relation © can be used 
to determine the free-energy difference. It is shown in [10) that the asymptotics 
(|12p generally satisfies the Crooks relation exactly, which we demonstrated for two 
representative examples of the potential (|22|) . 

Finally, we tested the universal form (|34|) of the tails of work distributions 
proposed by M. Palassini and F. Ritort against our results for the asymptotics (fl~2|) . 
For a broad range of parameters used for the model potential (f2"2"| . we found a 
good agreement between (|34|) and our asymptotics. Only if the work distribution 
differs markedly from a Gaussian form, a reliable fit of (|34j) is likely to require more 
data points than usually acquired in single-molecule experiments. This might lead 
to a significant difference between the exact and estimated free-energy difference as 
illustrated by our examples. 
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